Propagation of uncertainties in the Skyrme energy-density-functional model 
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Parameters of nuclear energy-density-functionals (EDFs) are always derived by an optimization 
to experimental data. For the minima of appropriately defined penalty functions, a statistical sen- 
sitivity analysis provides the uncertainties of the EDF parameters. To quantify theoretical errors 
of observables given by the model, we studied the propagation of uncertainties within the unedfO 
Skyrme-EDF approach. We found that typically the standard errors rapidly increase towards neu- 
tron rich nuclei. This can be linked to large uncertainties of the isovector coupling constants of the 
currently used EDFs. 
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Introduction. The nuclear cncrgy-density-functional 
(EDF) approach is the only microscopic method that can 
be applied throughout the entire nuclear chart. Owing 
to the existing and future radioactive-beam facilities, the 
experimentally known region of the nuclear landscape 
will be pushed towards more and more exotic systems. 
New experimental data on nuclei with a large neutron 
or proton excess will put predictions of the current and 
future EDF models to test. Therefore, it is important to 
assess the uncertainties and predictive power of models 

The nuclear EDF models can best be formulated in 
terms of effective theories, tailored to low-energy phe- 
nomena with high-energy effects taken into account 
through model parameter adjustments. Although a rig- 
orous application of the effective-theory methodology [5( 
is not yet fully implemented, the aspect of systematic ad- 
justments of model parameters to data is always a crucial 
element of the approach. Only very recently, such ad- 
justments are accompanied with the full covariance error 
analysis 

The aim of the present work is to study how uncer- 
tainties in the model parameters propagate to observ- 
ables calculated within the model. This appears to be 
the simplest and most fundamental method to assess the 
predictive power of the model. First, the model is char- 
acterized not only by the fit residuals, that is, deviations 
between measured and calculated observables, but also 
by propagated uncertainties, which carry the errors of 
model parameters over to calculated observables. In this 
way, even the data points that have been used for the pa- 
rameter adjustments are characterized by the propagated 
uncertainties. 

A properly balanced parameter adjustment should, in 
principle, lead to propagated uncertainties, which are 
comparable in magnitude to fit residuals. Second, ex- 
trapolations of the model to unmeasured data points, 
which is, in fact, one of the principal goals of formulat- 
ing models, bear no meaning if they are not accompanied 
by their propagated uncertainties. Indeed, large propa- 



gated uncertainties simply mean that relatively small de- 
viations of model parameters, within the bounds of the 
adjustment, strongly influence the predicted values. On 
the one hand, this may give us invaluable information on 
which parameters of the model are grossly undetermined, 
and on the other hand, this may indicate which observ- 
ables, when measured, would have the strongest impact 
on decreasing the errors of parameters. 

There is, of course, another important class of theoreti- 
cal uncertainties, which is not directly accessible through 
the statistical analysis of parameters, and which is very 
difficult to estimate, namely, the one stemming from the 
definition and limitation of the model itself, see discus- 
sion in Refs. 0, H, H|, ■ However, the statistical analysis 
may also here give us interesting information. Indeed, in 
cases when the propagated uncertainties turn out to be 
significantly smaller than the fit residuals, one may sus- 
pect that the model lacks important physics and/or new 
terms in the Hamiltonian (functional) and extended pa- 
rameter freedom. This is so, because we are then in the 
situation of model parameters well constrained by some 
data, but not describing some other data correctly. This 
can be contrasted with the situation of propagated uncer- 
tainties larger than the fit residuals, when one cannot tell 
whether the model is deficient or whether its parameters 
are underdetermined. 

Method. A good measure of theoretical uncertainty is 
to calculate the statistical standard deviation of an ob- 
servable, predicted by the model. In the present, work 
we used the Skyrme EDFs within the Hartree-Fock- 
Bogoliubov (HFB) nuclear structure calculations. To de- 
termine the propagation of uncertainties from the model 
parameters to some observable, predicted by the model, 
one needs information about the standard deviations and 
correlations between the model parameters. Some of the 
recent Skyrme-EDF parameter optimizations have pro- 
vided this information @,0,0- In this work, to determine 
the propagation of theoretical uncertainties, we have cho- 
sen to use the unedfO parameterization of the Skyrme 
EDF 0. 
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The calculated statistical error of an observable is 
linked with the covariance matrix of the nuclear EDF 
model parameters. From the Taylor expansion, we can 
find that the square of statistical error reads, 
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where y is the calculated observable, Xi is the parame- 
ter of the EDF and Cav(xi, Xj) is the covariance matrix 
between parameters x% and Xj. In the present work we 
use the covariance matrix that was obtained from the 
sensitivity analysis of the unedfO parameterization, see 
Tables VIII and X of Ref. Q ■ Another use of the covari- 
ance matrix consists in calculating covariance ellipsoids 
between pairs of observables [3 J8| and analyzing correla- 
tions between observables Qjl, [l3| . 

Calculation of the standard error ([I]) requires calcula- 
tion of derivatives of observable y with respect to model 
parameters Xj. As the relations between parameters 
and observables can be complicated, to calculate the 
derivatives one has to use numerical methods. By def- 
inition, the derivative of a function / can be approxi- 
mated simply by a finite difference, f'(a) m Go(h) = 
(f(a + h) — f(a — h)) /2h. To improve the accuracy, we 
used the Richardson extrapolation method [l4[ , 
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with m = 1, 2, . . ., which then satisfies, /'(a) — G m (h) = 
0(/j2(m+i)) Generally speaking, a larger m means im- 
proved accuracy on derivatives with an increased cost of 
computational time. In this work we extrapolate up to 
m = 2 as a balanced point. 

Large values h may lead to inaccurate derivatives, 
whereas small values may suffer from rounding errors. 
In this paper, we checked that when h changes for any 
parameter between 0.001% and 0.1% of the parameter's 
value, the calculated errors of observables do not change 
much. Thus values of h fixed in this region were used. 

The HFB calculations were carried out by using the 
computer code hosphe [ll| [lj| , which solves the HFB 
equations in a spherical harmonic-oscillator basis. We 
have run it within exactly the same conditions as those 
defined for the optimization of unedfO. In particular, 
calculations were done in the space of 20 major oscillator 
shells and the Lipkin-Nogami procedure of Ref. [13] with 
the pairing cut-off window of E cut = 60MeV was used. 
We determined ground-state properties of all even-even 
semi-magic nuclei with N = 20, 28, 50, 82, and 126 and 
Z = 20, 28, 50, and 82 extending between the two-proton 
and two-neutron driplines. 

For unedfO, some of the very neutron rich Ca, Sn, 
and Pb isotopes turn to be deformed, see, e.g., supple- 
mentary material of Ref. The same also holds for 
some semimagic isotonic chains at the very neutron-rich 
regime. Since the current calculations are done with a 
spherical HFB solver, the calculated uncertainties may 



lack deformation effects. A priori, one can expect that 
the standard errors determined in spherical and deformed 
nuclei should be similar, see also Ref. [lU, although in 
the future this conjecture must be tested in more exten- 
sive calculations. The aim of this work is to quantify a 
typical magnitude of the standard error related to an ob- 
servable, and thus assess the overall predictive power of 
the model. 

Results. In Fig. [TJ we present the calculated stan- 
dard errors of binding energies of semi-magic nuclei. The 
upper panel shows the isotonic chains with magic neu- 
tron numbers of N = 20, 28, 50, 82, and 126 and the 
lower panel shows the Ca, Ni, Sn, and Pb isotopic chains. 
Within the experimentally known regions of nuclei, where 
also data points were used during the EDF optimiza- 
tion, standard error usually stays rather constant and is 
typically of the order of 1 MeV. Moving to the neutron 
rich regime, calculated standard error increases rather 
steadily This is directly related to the fact that the 
isovector parameters of the EDF are not constrained as 
well as the isoscalar ones 0, Q • With isotonic (isotopic) 
chains, the nuclei with small Z (large N) of Fig. [T] are the 
ones with most neutron excess, and therefore standard 
deviation is largest in these nuclei. The values reached 
there can be of the order of 10 MeV, which means that 
for the total binding energies, when approaching the neu- 
tron drip line the predictive power of the model is rather 
poor. 

For the Pb and N = 126 nuclei we also show, as 
an illustration, the moduli of fit residuals with respect 
to experimental data [l!|. We note immediately, that 
the propagated error of the binding energy of 208 Pb 
of 0.41 MeV is much smaller than the fit residual of 
4.07 MeV. This clearly indicates that no refined fits can 
probably account for this experimental data point, be- 
cause the fit parameters are already quite rigidly con- 
strained by other data. This result points to a miss- 
ing physics in the description of the 208 Pb mass, which 
is most probably related to a poor description of the 
ground-states collective correlations in doubly magic sys- 
tems, see discussion in Refs. (20| - |22T |. 

In Fig. [21 the standard errors of two-neutron (/f>2n) and 
two-proton (S^p) separation energies are shown. Sim- 
ilarly as for the binding energies, within the experi- 
mentally know region, these standard errors stay rather 
small, that is, in heavy nuclei below 0.5 MeV for S2 P 
and below 0.25 MeV for S^n- Then, they start to in- 
crease rapidly when the calculations are extrapolated to 
neutron-rich regime, especially after the doubly magic 
gaps are crossed. Owing to the cancellation of errors, 
even in very exotic nuclei, the standard errors of two- 
particle separation energies are much smaller than those 
of the binding energies. This nicely illustrates the robust- 
ness of the position of two-neutron driplines discussed in 
Ref. [|. 

In Figs. [3] and SI we plot the standard errors of neutron 
and proton rms radii, respectively. One can see that the 
uncertainties related to the proton rms radii are much 



3 





N=20 




N=28 




N=50 


T 


N=82 




N=126 




Residue 





10 20 30 40 50 60 70 80 90 100 



10 20 30 40 50 60 70 80 90 100 




20 40 60 80 100 120 140 160 180 200 

N 

FIG. 1: (Color online) Standard errors of binding energies 
of semi-magic nuclei, calculated for (a) isotonic chains with 
magic neutron numbers and (b) isotopic chains with magic 
proton numbers. For the Pb and iV = 126 nuclei, moduli of 
fit residuals with respect to experimental data [lij are also 
shown. 



smaller than those corresponding to neutron rms radii. 
This is directly related to the fact that in the optimiza- 
tion of the unedfO parameter set, experimental input on 
the charge radii was used. For example, in 208 Pb we ob- 
tain cr(i? n )=0.062fm and cr(i? p )=0.006fm. The standard 
errors of neutron radii steadily increase towards the neu- 
tron drip line, again illustrating the large uncertainties 
in the isovector coupling constants of the functional. In 
heavy nuclei, the same is also true for protons, however, 
in light nuclei, the standard errors of proton radii also 
increase towards the proton drip line. This effect can 
be related to the weak binding of protons within a low 
Coulomb barrier. Uncertainties related to the proton and 
neutron rms radii are closely linked to the uncertainties 
of the neutron skin, see discussion in Ref. [IH . 

In Fig. [5j we present standard errors of the calculated 
mean-field (point-particle) proton and neutron densities 
in 208 Pb. Again we see that the standard errors of proton 
densities are significantly smaller than for neutrons, and 
that the neutron ones extend to much higher distances. 
The neutron density of Fig. [5ji can be compared with 
the experimental weak-charge profile presented in Fig. 1 
of Ref. [23|, obtained from the PREX measurement at 
Jefferson Lab (24[. In the region of relative flat neutron 
density profile (that is, R < 4fm) we find our standard 
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FIG. 2: (Color online) (a) Standard errors of two-proton sep- 
aration energies S^p for neutron semi-magic nuclei, (b) Stan- 
dard errors of two-neutron separation energies Sz n for proton 
semi-magic nuclei. 



errors to be 0.0035 fm~ or less. In the same region, the 
incoherent sum of experimental and model errors [23| is 
roughly 0.0075 fm~ 3 . Even though the density profile 
of Fig. [5] is for a point-neutron density, folding it with 
a finite weak-charge distribution of neutrons, or includ- 
ing in the experimental estimate other small corrections, 
should not significantly change our result of propagated 
theoretical uncertainties being about twice smaller than 
the experimental ones. 

Owing to the correlations between the model parame- 
ters, the standard errors calculated by using Eq. ^ can 
be strongly affected by the off-diagonal components of 
the covariance matrix. This is illustrated in Fig. [HI where 
various components of the sum contributing to squares of 
standard errors of the binding energy E (a), two- neutron 
separation energy S^n (b), neutron rms radius R n (c), 
and proton rms radius R p (d) in 208 Pb are plotted. The 
order of the parameters used in Fig. [6] is the same as 
that used in the definition of the correlation matrix of 
unedfO, see Tables VIII and X of Ref. Q. For the sake 
of completeness, these parameters are also listed in the 
caption of Fig. [5] 

One sees from Fig. [6^, that the standard error of the 
binding energy E results from a strong cancellation be- 
tween large positive (> 500 MeV 2 ) diagonal contributions 
coming from and (parameters No. 3 and 4), 

and large negative (< —500 MeV 2 ) contributions coming 
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FIG. 3: (Color online) Same as in Fig. [T] but for the neutron 
rms radii. 



FIG. 4: (Color online) Same as in Fig. [T] but for the proton 
rms radii. 



from the off-diagonal correlations between these two pa- 
rameters. The remaining parameters of the functional 
give much smaller positive and negative contributions, 
bringing the final value to cr 2 (£?) =0.1 7 MeV 2 . It is inter- 
esting to note that even in this well-bound nucleus, the 
largest uncertainties of the binding energy come from the 
poor symmetry-energy characterization of the functional. 

The situation obtained for the standard error of S2n 
(Fig. |6Jd) is quite similar. Although the value of 
CT 2 (S , 2n )=0.02MeV 2 is much smaller than that for the 
binding energy, it results form the cancellation of con- 
tributions coming form the above two parameters and 
C pAp (parameters No. 6), which is the parameter charac- 
terizing the isovector surface properties of the functional. 
Again, the uncertainties in the isovector characterization 
of the functional dominate the standard error here. 

The pattern obtained for the standard error of the 
neutron rms radius of cr 2 (i? n )=3.8x 10 -3 fm 2 (Fig. [Bt) is 
markedly different, with the single parameter dom- 
inating the uncertainty of this variable, although again 
the isovector properties are to be blamed. Only for the 
proton rms radius (Fig. |BJl), its fairly small standard er- 
ror of cr 2 (i?p)=4.1xl0 _5 fm 2 is given by contributions 
coming from and isoscalar parameters p c and Cq Ap . 

Even though the single-particle (s.p.) energies are 
not strictly-speaking observables, they are nevertheless 
closely connected to observable levels of neighboring 
even-odd nuclei and to other nuclear properties like, e.g., 
deformations of open-shell systems. It is, therefore, in- 
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FIG. 5: Standard errors of densities in 208 Pb for (a) protons 
and (b) neutrons as functions of the radial variable R. Panels 
(c) and (d) show the proton and neutron densities, respec- 
tively, with the standard errors represented by error bars. 



teresting to discuss their standard errors and check how 
tightly their values are fixed by the adjustment of the 
functional to other observables, see also discussion in 
Ref. 0. Since our calculations were performed with 
pairing correlations and Lipkin-Nogami corrections taken 
into account, we performed the analysis for the s.p. ener- 
gies e HF defined as the eigenstates of the mean-field part 
of the HFB Hamiltonian. 
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FIG. 6: Contributions of various components of the sum |T|) 
to squares of standard errors of the binding energy E (a) , two- 
neutron separation energy S21, (b), neutron rms radius R n (c), 
and proton rms radius R p (d) in 208 Pb. The numbers given in 
the ordinates and abscissas refer to the order of parameters 
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In Table|H we listed the s.p. energies and their standard 
errors calculated in 208 Pb. We see that the proton (neu- 
tron) rms of standard errors, equal to 0.196 (0.176) MeV, 
is significantly smaller than the rms deviations with re- 
spect to empirical data [H| of 0.52 (0.61) MeV. The rms 
deviations obtained in lead are still much smaller than 
those obtained when all doubly-magic nuclei are con- 
sidered simultaneously, which typically equal to about 
1.2 MeV (26[, quite independently of which functional 
and which data evaluation is used. Anyhow, we see again 
that the s.p. energies, which in the fit of the unedfO func- 
tional were not taken into account, are rather tightly con- 
strained by the fit to other observables. It is, therefore, 
unlikely that, within the parameter space of the standard 
Skyrmc functional, their agreement with data can be im- 
proved, which supports conclusions reached in Ref. [26j . 

Of course, in principle, the obtained insufficient agree- 
ment with data can also signal some important miss- 
ing physics. For the s.p. energies, an obvious and of- 
ten invoked effect is the well-known particlc-vibration- 
coupling (PVC) 27 1, which continues to be the subject 
of recent works 28l430i|. Before studying this problem 
for the Skyrme functionals in full detail [3l[ , in Table U 
we presented the PVC corrections to s.p. energies, 5e PVC , 
calculated with the standard methodology reviewed re- 
cently in, e.g., Ref. [29} . We see that in the particular 
case studied here, the inclusion of the PVCs does not 
improve the agreement with data, but, in fact, the rms 
deviation obtained for protons (neutrons) even increases 
to 0.70 (0.68) MeV. An improved agreement with data 
should probably be sought for by using richer parameter 
spaces of the functional, such as, e.g., those proposed in 
Refs. [jlHIl • Evidently, further studies of this problem 



TABLE I: The 208 Pb proton (top) and neutron (bottom) s.p. 
energies e H F (b) and their standard errors (c) , as compared to 
the empirical values (d) [2^], residuals, A(e H F) = e H F — e cxp 
(e), PVCs <5e PV c (f), and residuals of the PVC-corrected s.p. 



energies A(e PV c) = e H F + Se F 



e OX p (g). Where applica- 



ble, we also give the root-mean-square (rms) values of entries 
shown in a given column. All energies are in MeV. 



orbital 


e H F 


a(e HF ) 


6cxp 


A(e HF ) 


8e FVC 


A(epvc) 


(a) 


(b) 


(c) 


(d) 


(e) 


(f) 


(g) 


7T3pi/2 


0.219 


0.181 


-0.16 


0.38 


-0.440 


-0.06 


T3p 3 /2 


-1.045 


0.139 


-0.68 


-0.37 


-0.662 


-1.03 


7r2f 5 /2 


-1.284 


0.229 


-0.97 


-0.31 


-0.480 


-0.79 


Tlil3/2 


—2.794 


0.238 


-2.10 


-0.69 


-0.221 


-0.92 


7rlh 9/ 2 


-3.501 


0.282 


—3.80 


0.30 


-0.280 


0.02 


7r2f 7 / 2 


-3.725 


0.115 


-2.90 


-0.83 


-0.284 


—1.11 


7t3si /2 


-8.036 


0.140 


—8.01 


—0.03 


—0.108 


-0.13 


7r2d 3 /2 


-8.378 


0.199 


-8.36 


-0.02 


0.220 


0.20 


Tlhn/2 


-9.153 


0.207 


-9.36 


0.21 


-0.141 


0.07 


7r2d 5/ 2 


—10.117 


0.117 


—9.82 


—0.30 


0.116 


-0.18 


flgr/2 


-10.908 


0.229 


-12.00 


1.09 


0.131 


1.22 


rms 


n.a. 


0.196 


n.a. 


0.52 


0.328 


0.70 


^3d 3/2 


-1.856 


0.250 


-1.40 


-0.46 


-0.104 


-0.56 


u4s 1/2 


-2.051 


0.235 


-1.90 


-0.15 


-0.668 


-0.82 


v2g 7/2 


-2.141 


0.258 


-1.44 


-0.70 


-0.280 


-0.98 


T/lfi K /o 

" - Li lo/ A 


-2.231 


0.167 


-2.51 


0.28 


-0.226 


0.05 


!/3d 5/ 2 


-2.549 


0.166 


-2.37 


-0.18 


-0.384 


-0.56 


^lin/2 


-2.680 


0.223 


-3.16 


0.48 


-0.271 


0.21 


^2g 9/2 


-4.336 


0.071 


-3.94 


-0.40 


-0.183 


-0.58 


i/3pi /2 


-7.855 


0.109 


-7.37 


-0.49 


0.152 


-0.33 


^3p 3/2 


-8.503 


0.065 


-8.26 


-0.24 


-0.119 


-0.36 


^2f 5/2 


-8.519 


0.143 


-7.94 


-0.58 


0.156 


-0.42 


^lil3/2 


-8.603 


0.162 


-9.24 


0.64 


-0.130 


0.51 


i/lh 9 /2 


-9.922 


0.190 


-11.40 


1.48 


0.120 


1.60 


^2f 7/2 


-10.407 


0.103 


-9.81 


-0.60 


0.179 


-0.42 


rms 


n.a. 


0.176 


n.a. 


0.61 


0.273 


0.68 



are very much required. 

Summary and Conclusions. In the framework of the 
unedfO Skyrme-EDF model, we studied propagation of 
uncertainties from model parameters to calculated ob- 
servables. We found that the magnitude of standard 
errors increases towards neutron rich nuclei. This is di- 
rectly linked to the isovector part of the functional, where 
coupling constants have large uncertainties. 

For binding energies, the standard errors are in the 
range of one to several MeV - towards the neutron drip 
line rapidly increasing to about 10 Mev. For two-particle 
separation energies, they range from a few hundred keV 
to a few MeV. Compared to the deviations from experi- 
mental values, the standard errors of binding energies are, 
in the experimentally known region, of the same order of 
magnitude or somewhat smaller. For two-particle sepa- 
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ration energy, standard errors are typically smaller, but 
still comparable. It is worth stressing that the theoretical 
uncertainties of binding energies and two-particle sepa- 
ration energies are still much larger than the present ex- 
perimental errors in Penning-trap measurements [33l |34| . 
For proton rms radii, our standard errors are also compa- 
rable, but slightly smaller than the rms deviations from 
experimental values. Theoretical uncertainties related to 
the neutron skin are actually smaller than the current 
experimental estimates [T§ |. 

In conclusion, determination of theoretical uncertain- 
ties in terms of standard statistical errors of observables 
provides a way to assess the predictive power of nuclear 
models. In the future EDF parameter optimizations, 
it would be useful to perform by default the sensitivity 



analyses of the obtained minima. With this information, 
theoretical uncertainties related to the model parameters 
and predictions should be routinely assessed. 
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